IJGI

First, we load the data into the workspace.

Following is a summary of the dataset.

summary(checkin.NY)

##    weather.id        gid             user_id            gender      
##  Min.   :  24   Min.   :      1   Min.   :      32   female:259667  
##  1st Qu.: 888   1st Qu.: 276193   1st Qu.: 1542768   male  :311595  
##  Median :1763   Median : 552188   Median : 6913375   NA's  :  8524  
##  Mean   :1760   Mean   : 551507   Mean   :15589057                  
##  3rd Qu.:2589   3rd Qu.: 825534   3rd Qu.:21767105                  
##  Max.   :3590   Max.   :1135513   Max.   :89322077                  
##                                                                     
##                      venue_id           lat            lon       
##  43a52546f964a520532c1fe3: 12663   Min.   :40.5   Min.   :-74.2  
##  4ace6c89f964a52078d020e3:  7886   1st Qu.:40.7   1st Qu.:-74.0  
##  42911d00f964a520f5231fe3:  3521   Median :40.7   Median :-74.0  
##  4ae6363ef964a520aba521e3:  2967   Mean   :40.7   Mean   :-74.0  
##  49b7ed6df964a52030531fe3:  2835   3rd Qu.:40.8   3rd Qu.:-74.0  
##  42829c80f964a5206a221fe3:  2191   Max.   :40.9   Max.   :-73.7  
##  (Other)                 :547723                                 
##   timestamps.x                             cate_l1      
##  Min.   :1.39e+09   Food                       :170208  
##  1st Qu.:1.39e+09   Shop & Service             : 90591  
##  Median :1.40e+09   Nightlife Spot             : 70389  
##  Mean   :1.40e+09   Travel & Transport         : 62418  
##  3rd Qu.:1.40e+09   Arts & Entertainment       : 60727  
##  Max.   :1.40e+09   Professional & Other Places: 57076  
##                     (Other)                    : 68377  
##                  cate_l2          datetime                  
##  Gym / Fitness Center: 32386   Min.   :2014-02-01 00:00:04  
##  Bar                 : 25039   1st Qu.:2014-03-09 19:27:56  
##  Airport             : 24891   Median :2014-04-15 11:12:26  
##  Office              : 24356   Mean   :2014-04-15 06:30:21  
##  Coffee Shop         : 18346   3rd Qu.:2014-05-20 03:09:42  
##  American Restaurant : 15471   Max.   :2014-06-30 23:58:52  
##  (Other)             :439297   NA's   :69                   
##       hour          yearday            weekday            isweekend     
##  19     : 57300   Length:579786      Length:579786      Weekend:192831  
##  18     : 52049   Class :character   Class :character   Workday:386955  
##  20     : 45170   Mode  :character   Mode  :character                   
##  13     : 38504                                                         
##  17     : 37898                                                         
##  14     : 35523                                                         
##  (Other):313342                                                         
##     last.gid       last.timestamps                      last.venue_id   
##  Min.   :     -1   Min.   :-1.00e+00   43a52546f964a520532c1fe3:  9083  
##  1st Qu.: 203258   1st Qu.: 1.39e+09   4ace6c89f964a52078d020e3:  6081  
##  Median : 491981   Median : 1.40e+09   42911d00f964a520f5231fe3:  2972  
##  Mean   : 498473   Mean   : 1.29e+09   4ae6363ef964a520aba521e3:  2628  
##  3rd Qu.: 784803   3rd Qu.: 1.40e+09   49b7ed6df964a52030531fe3:  2456  
##  Max.   :1135513   Max.   : 1.40e+09   (Other)                 :513535  
##                                        NA's                    : 43031  
##     last.lat       last.lon                          last.cate_l1   
##  Min.   :-1.0   Min.   :-74.2   Food                       :158896  
##  1st Qu.:40.7   1st Qu.:-74.0   Shop & Service             : 86910  
##  Median :40.7   Median :-74.0   Nightlife Spot             : 65845  
##  Mean   :37.6   Mean   :-68.5   Arts & Entertainment       : 54447  
##  3rd Qu.:40.8   3rd Qu.:-73.9   Professional & Other Places: 53443  
##  Max.   :40.9   Max.   : -1.0   (Other)                    :117214  
##                                 NA's                       : 43031  
##                last.cate_l2     last.user_id      time.interval   
##  Gym / Fitness Center: 31436   Min.   :      32   Min.   :     0  
##  Bar                 : 23525   1st Qu.: 1541303   1st Qu.:     3  
##  Office              : 22786   Median : 6812580   Median :    21  
##  Airport             : 18457   Mean   :15394230   Mean   : 28852  
##  Coffee Shop         : 17437   3rd Qu.:21345141   3rd Qu.:    93  
##  (Other)             :423114   Max.   :89322077   Max.   :390052  
##  NA's                : 43031   NA's   :43031                      
##     latitude      longitude                  conds           winddird    
##  Min.   :40.7   Min.   :-73.9   Clear           :289609   Min.   :  0.0  
##  1st Qu.:40.7   1st Qu.:-73.9   Overcast        :134203   1st Qu.:  0.0  
##  Median :40.7   Median :-73.9   Mostly Cloudy   : 46510   Median :  0.0  
##  Mean   :40.7   Mean   :-73.9   Partly Cloudy   : 31744   Mean   : 92.7  
##  3rd Qu.:40.7   3rd Qu.:-73.9   Scattered Clouds: 28311   3rd Qu.:190.0  
##  Max.   :40.7   Max.   :-73.9   (Other)         : 47843   Max.   :360.0  
##                                 NA's            :  1566                  
##     windspd        temperatur       fog             rain        
##  Min.   : 0      Min.   :-12.2   Mode :logical   Mode :logical  
##  1st Qu.: 6      1st Qu.:  4.4   FALSE:576316    FALSE:551962   
##  Median : 9      Median : 12.8   TRUE :3470      TRUE :27824    
##  Mean   :10      Mean   : 11.8   NA's :0         NA's :0        
##  3rd Qu.:13      3rd Qu.: 19.4                                  
##  Max.   :35      Max.   : 31.1                                  
##  NA's   :38929   NA's   :188                                    
##     snow          thunder         tornado       
##  Mode :logical   Mode :logical   Mode :logical  
##  FALSE:573449    FALSE:579786    FALSE:579786   
##  TRUE :6337      NA's :0         NA's :0        
##  NA's :0                                        
##                                                 
##                                                 
##

1. Exploratory analysis

Temporal: Frequency domain

freq.plot(checkin.NY,title="New York City")

plot of chunk unnamed-chunk-4

The analysis in the frequency domain reveals a strong harmonics with a period of 24 hours. Therefore, it is reasonable to analyze the data by 24 hours.

Temproal: Time domain

The distribution of the check-in categories accross one period (24 hours).

time.distribution.plot(checkin.NY,title="New York City")

plot of chunk unnamed-chunk-5

Temporal: Radial plot

The distribution of the check-in categories in 24 hours in radial plots.

time.radial.plot(checkin.NY)

plot of chunk unnamed-chunk-6

Spatial

have a look at a intuitive spatiotemporal visualizaiton of the data…

basemap = map.plot("../../global/data/shapefiles", 
                   "NYC_borough_boundaries_WGS84",
                   alpha=0.1,size=0.3,color="grey")

saveGIF({
  point.animation.plot(checkin.NY,title="New York City", 
                       color="cate_l1", color.concrete=FALSE,
                       basemap=basemap,axis.size=18,
                       title.size=24,legend.size=18,point.size=3)
}, interval = 0.5, movie.name = "spatiotemporal.gif", 
ani.width = 1500, ani.height = 1200)

What will the statistics say when we neglect the temporal factors? First of all, if we try to find spatial clusters based on differnt checkin categories, the distribution of the founded clusters are quite differnt. It indicates that each category has differnt correlations with the geographic space.

if(!exists("basemap")){
    basemap <- map.plot("../../global/data/shapefiles", 
                   "NYC_borough_boundaries_WGS84",
                   alpha=0.1,size=0.3,color="grey")
}

## OGR data source with driver: ESRI Shapefile 
## Source: "../../global/data/shapefiles", layer: "NYC_borough_boundaries_WGS84"
## with 5 features and 4 fields
## Feature type: wkbMultiPolygon with 2 dimensions

plots <- lapply(split(checkin.NY,checkin.NY$cate_l1),function(ci){
    # find out clusters for the type of category
    clusters = spatial.clustering(ci)
    
    centers = clusters[["centers"]]
    points = clusters[["point.unique"]]
    wss = clusters[["wss"]]
    pct = clusters[["pct"]]
    lbl = paste(nrow(centers),"clusters;\nWSS:",
                  formatC(wss,digits=2,format="f"),
                  "(",format.percent(pct),")")
    
    # add basic points
    gg.map <- point.plot(points,x="lon.x",y="lat.x",alpha=0.05, basemap=basemap)
    
    # add cluster information
    gg.map <- point.plot(centers, x="lon.center",y="lat.center",
                         color= "cid.ordered", color.concrete=FALSE,
                         point.size = log(centers$size,5),alpha = 0.7,
                         basemap=gg.map,
                         xlim=range(checkin.NY$lon,finite=TRUE),
                         ylim=range(checkin.NY$lat,finite=TRUE))
    
    # some plot configuration
    gg.map <- gg.map + ggtitle(ci[1,"cate_l1"]) +
        theme(legend.position="none") + 
        geom_text(aes(x = -74.13, y = 40.87), label = lbl, size=2) 
    
})

## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt

# png("categorized_clustering.png",width=15*ppi,height=6*ppi,res=ppi)
# grid.arrange(plots[[1]],plots[[2]],plots[[3]],plots[[4]],
#              plots[[5]],plots[[6]],plots[[7]],plots[[8]], 
#              plots[[9]],plots[[10]],nrow=4, ncol=3)
# dev.off()

plots[[1]];plots[[2]];plots[[3]];plots[[4]];plots[[5]];plots[[6]];plots[[7]];plots[[8]]; plots[[9]];plots[[10]]

plot of chunk unnamed-chunk-8

We could also see the most domimant categories in each part of the entire city.

# prepare polygon data
SPDF = readOGR(dsn = "../../global/data/shapefiles", layer = "NYC_zipcode")

## OGR data source with driver: ESRI Shapefile 
## Source: "../../global/data/shapefiles", layer: "NYC_zipcode"
## with 236 features and 6 fields
## Feature type: wkbPolygon with 2 dimensions

# intersect the checkin data with the each postal code region
checkin.NY = point.in.poly(checkin.NY, SPDF, copy.attr="POSTAL")

## [1] "Spatial data has been prepared."
## [1] "The information from the polygon has been assigned to the points."

# get the distribution of categories in each postal code region
cate.distr=cate.distr.in.poly(checkin.NY)

## [1] "The distribution of category has been assigned to each polygon."

# and use that information to implement the SPDF
SPDF = poly.with.category(SPDF, cate.distr=cate.distr, poly.attr="POSTAL")

# plot

mapdf=df.from.spdf(SPDF)
# mapdf$density=apply(mapdf,1,function(i){i[i["cate.dom"]]})
# mapdf$density=as.numeric(formatC(mapdf$density,digits=1,format = "f"))
# png("main_category.png",height=6*ppi, width=16*ppi,res=ppi)
# grid.arrange(
map.plot(mapdf=mapdf,
         more.aes=aes_string(fill="Category.1st"),
         color="grey",size=0.3,alpha=0.9)+
    theme_bw()+
    xlab("")+ylab("")+
    theme(legend.title = element_text(size=8),
          legend.text = element_text(size = 8))

## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt

plot of chunk unnamed-chunk-9

map.plot(mapdf=mapdf,
         more.aes=aes_string(fill="Category.2nd"),
         color="grey",size=0.3,alpha=0.7)+
    theme_bw()+
    xlab("")+ylab("")+
    theme(legend.title = element_text(size=8),
          legend.text = element_text(size = 8))

## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt

plot of chunk unnamed-chunk-9

# nrow=1, ncol=2, widths=c(1,1) )
# dev.off()

Meteorological: Correspondence analysis

cate.conds = xtabs(~conds+cate_l1, data=checkin.NY)
#prop.table(cate.conds, 1) # row percentages
#prop.table(cate.conds, 2) # column percentages
fit <- ca(cate.conds)
#print(fit) # basic results
summary(fit) # extended results

## 
## Principal inertias (eigenvalues):
## 
##  dim    value      %   cum%   scree plot               
##  1      0.003244  52.2  52.2  *************************
##  2      0.002082  33.5  85.8  ****************         
##  3      0.000511   8.2  94.0  ****                     
##  4      0.000147   2.4  96.4  *                        
##  5      7.6e-050   1.2  97.6  *                        
##  6      6.7e-050   1.1  98.7                           
##  7      5.1e-050   0.8  99.5                           
##  8      2.4e-050   0.4  99.9                           
##  9      8e-06000   0.1 100.0                           
##         -------- -----                                 
##  Total: 0.006210 100.0                                 
## 
## 
## Rows:
##      name   mass  qlt  inr    k=1 cor ctr    k=2 cor ctr  
## 1  | Cler |  501  779   88 |   22 452  76 |  -19 327  86 |
## 2  |  Fog |    2  634   31 | -265 628  37 |  -26   6   1 |
## 3  | Haze |   22  975  178 | -200 798 272 |   94 177  94 |
## 4  | HvyR |    5  917   75 |  208 421  61 |  226 496 111 |
## 5  | HvyS |    0  116    5 | -170  20   0 |  377  96   1 |
## 6  | LgFR |    0  714   29 | -590 712  40 |   28   2   0 |
## 7  | LghR |   33  947   88 |   77 362  61 |   98 586 154 |
## 8  | LghS |    7  246   21 |  -62 203   8 |   29  44   3 |
## 9  | MstC |   80  738   38 |  -31 319  23 |  -35 420  47 |
## 10 | Ovrc |  232  957  107 |  -25 217  45 |   46 740 237 |
## 11 | PrtC |   55  592   41 |   38 315  25 |  -36 277  34 |
## 12 | Rain |   10  943   67 |  163 632  81 |  114 312  62 |
## 13 | SctC |   49  874   97 |  -59 287  53 |  -85 587 169 |
## 14 | Snow |    4  850  134 | -421 849 218 |    4   0   0 |
## 
## Columns:
##      name   mass  qlt  inr    k=1 cor ctr    k=2 cor ctr  
## 1  | ArtE |  105  763   76 |   49 526  76 |   33 237  54 |
## 2  | CllU |   17  677   27 |  -75 580  30 |   31  97   8 |
## 3  | Evnt |    5  557   20 |  -24  23   1 | -113 534  32 |
## 4  | Food |  294  688   48 |   18 301  28 |  -20 388  56 |
## 5  | NghS |  122  955  206 |   93 826 325 |   37 129  80 |
## 6  | OtdR |   69  890  200 |   -9   4   2 | -127 886 528 |
## 7  | PrOP |   98  986  272 | -123 878 457 |   43 108  87 |
## 8  | Rsdn |   27  352   42 |   23  54   4 |   54 298  38 |
## 9  | ShpS |  156  661   49 |  -29 430  41 |  -21 230  34 |
## 10 | TrvT |  107  791   59 |  -33 318  36 |   40 473  84 |

#plot(fit) # symmetric map
plot(fit, mass = TRUE, contrib = "absolute", map ="rowgreen", 
     arrows = c(TRUE, FALSE)) # asymmetric map

plot of chunk unnamed-chunk-10

3. Model - derivation and corresponding functions

Under the assumption \(H=i\) is independent from \(W=j\) \[P(C=k|H=i,W=j)=\frac{P(C=k,H=i,W=j)}{P(H=i,W=j)}=\frac{P(H=i,W=j|C=k)*P(C=k)}{P(H=i)*P(W=j)} (1) \]

since \(H=i\) is independent from \(W=j\), \[Exp[P(H=i,W=j|C=k)]=P(H=i|C=k)*P(W=j|C=k) (2)\]

therefore, \[Exp[P(C=k|H=i,W=j)]=Exp[ \frac{P(H=i,W=j|C=k)*P(C=k)} {P(H=i)*P(W=j)}] \\\ =\frac{P(H=i|C=k)*P(W=j|C=k)*P(C=k)}{P(H=i)*P(W=j)} \\\ =\frac{\frac{P(H=i,C=k)}{P(C=k)}*\frac{P(W=j,C=k)}{P(C=k)}*P(C=k)}{P(H=i)*P(W=j)} \\\ =\frac{P(C=k|H=i)*P(H=i)*P(C=k|W=j)*P(W=j)}{P(H=i)*P(W=j)*P(C=k)} \\\ =\frac{P(C=k|H=i)*P(C=k|W=j)}{P(C=k)} (3)\]

relevance contextualized by temporal factor

\[P_{u}(C=k|H=i)=\frac{\Phi_{u}(C=k,H=i)}{\Phi_{u} (H=i)} (4)\]

get.temporal.impact <- function(dataframe,hour){
    # dataframe.in.hour = checkin.single[which(checkin.single$hour==hour),]
    dataframe.in.hour = dataframe[which(dataframe$hour==hour),]
    phi.h = nrow(dataframe.in.hour)
    
    list.category = split(dataframe.in.hour, dataframe.in.hour$cate_l2)
    sapply(list.category, function(i){
        nrow(i)/phi.h
    })
}

relevance contextualized by (unweighted meteorological) factor

\[P_{u}(C=k|W=j)=\frac{Intercept(C=k,W=j)}{\sum Intercept(C,W=j)} (5)\]

get.meteorological.impact <- function(fit,conds){
    
    conds.id = which(fit[["rownames"]]==conds)
    ref.vec = fit[["rowcoord"]][conds.id,1:8]
    cate.all = fit[["colcoord"]][,1:8]
    
    intercepts = apply(cate.all, 1, function(x){ 
        (x[1]*ref.vec[1] + x[2]*ref.vec[2] + x[3]*ref.vec[3] + x[4]*ref.vec[4] +
             x[5]*ref.vec[5] + x[6]*ref.vec[6] + x[7]*ref.vec[7] + x[8]*ref.vec[8] ) /
            (ref.vec[1]^2 + ref.vec[2]^2 + ref.vec[3]^2 + ref.vec[4]^2 +
                 ref.vec[5]^2 + ref.vec[6]^2 + ref.vec[7]^2 + ref.vec[8]^2 ) 
        } )
    
#     vec = intercepts / sum(intercepts)  !!!wrong!!! intercetp can be negative
    vec = (intercepts - min(intercepts)) / (max(intercepts)-min(intercepts)) # scale into [0,1]
    names(vec) = fit[["colnames"]]
    
    vec
}

get.meteorological.impact2 <- function(dataframe,conds){
    
    dataframe.in.conds = dataframe[which(dataframe$conds==conds),]
    phi.c = nrow(dataframe.in.conds)
    
    list.category = split(dataframe.in.conds, dataframe.in.conds$cate_l2)
    sapply(list.category, function(i){
        nrow(i)/phi.c
    })
}

weighted meteorological factor

\[P_{u}^{*} (C=k|W=j)= w_{j}*[P_{u}(C=k|W=j)-\bar P_{u}]+\bar P_{u}\]

get.weather.weight <- function(fit){
    
    conds.all = fit[["rowcoord"]][,1:2]
    
    mag = apply(conds.all, 1, function(x){ 
        sqrt( (x[1]^2+x[2]^2) )
        } )
    
    mag / sum(mag)
}

get.weighted.meteorological.impact <- function(fit, conds, weight){
    
#     cate.conds = xtabs(~conds+cate_l2, data=dataframe)
#     fit <- ca(cate.conds)
    
    unweighted = get.meteorological.impact(fit,conds)
#     weights = get.weather.weight(fit)
    
#     conds.id = which(fit[["rownames"]]==conds)
    
#     vec = weights[conds.id] * (unweighted-mean(unweighted)) + mean(unweighted)
    vec = weight * (unweighted-mean(unweighted)) + mean(unweighted)
    names(vec) = fit[["colnames"]]
    
    vec
    
}


get.weighted.meteorological.impact2 <- function(dataframe, conds, weight){
    
    
    unweighted = get.meteorological.impact2(dataframe, conds)

    weight * (unweighted-mean(unweighted)) + mean(unweighted)
    
    
}

denominator (unweighted relevance)

\[P(C=k)=\frac{\Phi_{u} (C=k) }{\Phi_{u} }\]

get.denominator <- function(dataframe){
    
    phi.h = nrow(dataframe)
    
    list.category = split(dataframe, dataframe$cate_l2)
    denominator = sapply(list.category, function(i){
        nrow(i)/phi.h
    })
    
    denominator

    
}

the final result

\[E[P(C=k|H=i,W=j)]=\frac{P(C=k|H=i)*P(C=k|W=j)}{P(C=k)}\]

update @2014.09.04

This equation has been somehow verified :
- although p.k has a relatively high performance, multipling it with others do not bringer a higher one; instead, dividing does.
There should be some problems in the meteorological impacts.
- The performance based on calculation of interception is so poor that it totally brings down the whole performance. Instead, if just the simple probability distribution is adopted, the performance is getting better.

update @2014.09.05

problem has been found out. In the earlier function get.meteorological.impact(), I used intercept/sum(intercept), where intercept can be negative, and the sum(intercept) can be nonsense. Now it is fixed to (intercepts - min(intercepts)) / (max(intercepts)-min(intercepts)) to scale into [0,1].
besides, about the equation again:
- if p.kj is calculated by just using the probability distribution, then the equation p.kij = p.ki * p.kj / p.k is correct (which confirmes our derivation);
- if p.kj is calculated by personal CA, it should be noticed that the scaled interception does not describe the probability of C=k under the condition W=j. Instead, it describes the negative/positive impacts. Hence, the conditional probability should be the scaled interception multipled by the priori p.k. Therefore, here the final result should be p.kij = p.ki * p.kj.
- with a little differnt equations and different methods of getting p.j, the above two ways have very similar results (which is expected!).